Every timestep Δt, assumed equal to 1 hour, the power balance between the pre-heating power demand of natural gas W_gas, for the given input conditions and output set points, and the heat output that can be supplied by the heat pump W_HP,th is calculated.

The heat pump always has priority whenever there is a surplus of renewable thermal energy (i.e., electric power from the PV), all the pre-heating requirements are fulfilled with the heat pump, and the equivalent surplus electricity is sold to the grid. On the other hand, if the heat output supplying the HP is zero or insufficient, the auxiliary boiler comes into action, and the necessary NG flow rate is obtained from the primary flow via the Splitting Valve (SPV).

The control balance of the logic is described by equation (1) where W_BU is the boiler power to cover the thermal power deficit W_def, W_grid is the electrical power fed into the grid by the photovoltaic surplus W_surp; the various models of the components of the equation will be described in the following paragraphs.

$\Delta W(t)=W_{gas}(t)-W_{th,HP}(t)=\left\{\begin{array}{l}0\to W_{BU}=0;W_{grid}=0\\ <0\to W_{BU}=0;W_{grid}=W_{surp}\\ >0\to W_{BU}=W_{def};W_{grid}=0\end{array}\right.$ (1)

The electric power which can be exploited by the ASHP every hour W_el,HP is obtained following a control logic that compares the power output of the solar panels W_PV(t) with a “cut-off” threshold ε: when the solar field output is equal or higher to the ε the heat pump is switched on, but when the value of W_el,PV(t) falls below ε the heat pump is switched off.

$W_{el,HP}(t)=\left\{\begin{array}{l}{W}_{el,PV}(t)if{W}_{el,PV}(t)\ge \epsilon \\ 0if{W}_{el,PV}(t)<\epsilon \end{array}\right.$ (2)

The cut-off threshold is obtained by considering the advice given by the DSO and is set to 1 kW to avoid drops in the heat pump’s efficiency below the value of supplied electrical power. The heat output that the heat pump can provide every hour will be given by the product of the available electrical power and the actual coefficient of performance COP(t).

$W_{th,HP}(t)=COP(t)\times W_{el,HP}(t)$ (3)

Every timestep Δt, if the total power balance is higher than 0 (thermal energy deficit), the heat output to be supplied by the auxiliary boiler is calculated according to (4) as the difference between the required heat output and the heat output supplied by the HP calculated with the previous equation.

$W_{BU}(t)=W_{gas}(t)-W_{th,HP}(t)$ (4)

The total annual thermal energy supplied by the auxiliary boiler [kWh/year] is given by:

$E_{gas,y}={\displaystyle \sum _{i=1}^{N_{days}}{W}_{BU}(t)\times \Delta t}$ (5)

The total annual thermal energy saved will be the thermal energy supplied by the pump instead of, or together with, the auxiliary boiler.

$E_{th,sav,y}={\displaystyle \sum _{i=1}^{N_{days}}{W}_{th,HP}(t)\times \Delta t}$ (6)

When it comes to RES systems based on PV plants, the main parameters to be considered to assess the self-sufficiency level of the system are the SSR (Self Sufficiency Ratio) and SCR (Self Consumption Ratio), generally defined as ratios of amounts of electricity [17].

In this study, since annual demand is thermal, the SSR of the RES system is adjusted to the considered case study and defined as the ratio between the self-consumed thermal energy and the total yearly energy demand (7); these two parameters are obtained from the previous equations.

$SSR=\frac{E_{th,sav,y}}{E_{gas,y}}$ (7)

On the other hand, the SCR can be expressed as the ratio of self-consumed electric energy and the total yearly energy production.

$SCR=\frac{E_{el,sav,y}}{E_{PV,y}}=\frac{E_{el,sav,y}}{E_{el,sav,y}+E_{grid,y}}$ (8)

Where E_el,sav,y is obtained directly from the E_th,sav,y knowing the actual COP value every timestep according to environmental conditions, E_PV,y is the annual PV electric energy output and E_grid,y is the annual electricity sold to the grid.

The temperature of the gas arriving at the CGS is generally unknown because of the lack of sensors at the station. Therefore, it is assumed constant in several scientific works and equal to the worst case possible [11], [12], i.e. 0 °C. In this work, a more realistic model is exploited [8], which calculates the soil temperature surrounding a pipe buried at a 1-metre depth close to the CGS according to the variation of the air ambient temperature, assuming the NG temperature inside the pipes equal to the soil temperature (18).

$T_{soil}=0.0084\times T_{env}^2+0.3182\times T_{env}+11.403$ (17)

$T_{gas,in}=T_{soil}$ (18)

Figure 2 shows the NG inlet temperature increasing with the outside air temperature according to the model given by equation (17). For outside air temperature variations between 0 °C and 20 °C, the gas temperature varies between 11 °C and 21 °C. As the outside temperature decreases and drops below zero, the ground temperature stays at a minimum of around 10 °C.

Figure 2.

Natural gas inlet temperature vs. air ambient temperature according to (17)

The temperature change the gas undergoes during an adiabatic expansion depends on the final and initial pressure states and how the expansion is carried out. In a free expansion, the gas does no work and absorbs no heat, so the internal energy is conserved. Expanding in this way, the temperature of an ideal gas would remain constant, but the temperature of a real gas decreases, except at very high temperatures. On the other hand, the Joule-Thomson expansion method is intrinsically irreversible. During this expansion, the enthalpy remains unchanged, but unlike a free expansion, work is done that causes a variation in internal energy. This change due to the irreversibility of the process means that much greater cooling or heating can be achieved than in the case of free expansion.

The Joule-Thomson effect is a phenomenon whereby the temperature of a real gas decreases following expansion conducted at constant enthalpy. In literature, this effect is often assumed constant and equal to 4−5 °C/MPa during the gas throttling process inside a CGS [11], [12], or its calculation is avoided by imposing an isenthalpic transformation between the starting point and the end point [13]. For this work purpose, the coefficient is calculated according to the following formula [19], [20]:

${\mu }_{JT}={\left(\frac{\partial T}{\partial P}\right)}_H$ (19)

The Joule-Thomson coefficient μ_JT is calculated for several NG mixtures stored in the Cool Prop database [21], which are enlisted and described in Table 1. The gas outlet condition (P_out, T_out) is kept fixed and equal to the ideal set point values for pressure and several set point temperatures (30 MPa, 10 °C) to compute the isenthalpic process necessary for the μ_JT evaluation.

Figure 3 shows the linear dependence between the Joule-Thomson coefficient, which varies approximately between 4.1 °C/MPa and 5.2 °C/MPa, and the inlet pressure for all the considered NG mixtures. Pure methane (100% CH4) is fluid with the lowest value of μ, while the Ekofisk (North European) NG is the mixture with the highest temperature drop during the throttling phase at constant enthalpy.

Figure 3.

Joule-Thomson coefficient μ_JT for several gas mixtures vs. the inlet pressure P_in

The typical NG μ_JT curve is chosen to be the one that will be used in the following chapters to consider a general NG mixture composition since, in Italy, there is a very high variability of gas composition due to the multitude of import origins. Focusing on the typical NG μ-curve, Figure 4 highlights the effect of a different output temperature set point on the μ_JT coefficient and thus on the final pre-heating energy demand. As the set point at the gas outlet increases, the temperature drops to be made by the gas at the same inlet pressure at the CGS are reduced, and thus the difference between the inlet temperature and the outlet temperature ΔT_OI calculated with (11). On the other hand, a higher set point leads to an increase in the pre-heating factor calculated with (12).

Figure 4.

Joule-Thomson coefficient μ_JT for the natural gas typical composition for several set points of gas outlet temperature vs. the inlet pressure P_in

In this article, a real dataset covering one year of operation of a CGS in central Tuscany will be exploited. Figure 5a and Figure 6 show the distribution curves for the gas flow rate, the arrival pressures from the transport network, the external air temperature, and the gas outlet temperature.

Figure 5.

Gas volumetric flow rate (a) and inlet pressure (b) examples for the considered case study

Figure 6.

Gas temperature distribution at valve inlet (a) and valve outlet (b)

The lines of the graphs were obtained by plotting kernel density estimates to smooth the distribution and show the trend of the univariate variable. Figure 5 shows that the variation of the two input quantities, i.e., flow rate and pressure, is very wide. It is, therefore, important to consider this variation by giving the correct flow and pressure inputs to the model in equation (1).

Figure 6, on the other hand, shows the distributions of ambient temperature and, thus, of gas inlet temperature at the CGS and outlet temperature. The variation in ambient temperature is very wide and ranges from a few degrees below zero to temperatures of over 35 °C. It is interesting to assess how the gas outlet temperature follows the set point temperature set very precisely, slightly overheating compared to the set point.

The gas outlet temperature T_gas,out is the key parameter to be monitored, and it is strictly dependent on the value of the outlet temperature set point, which the DSO sets inside all the CGS. According to the Italian framework, this value must equal or exceed 5 °C. Still, for safety reasons, it is generally set at higher values and can be modulated for two working seasons: winter and summer. In this work, two different values for the outlet gas temperature are considered, replicating the actual output temperature setting inside a CGS in Italy, as can be seen in the following equation:

$T_{gas,out}=T_{gas,SP}=\left\{\begin{array}{l}{8}^{\circ }\text{C}\text{(}\mathit{\text{winter}}season)\\ {10}^{\circ }\text{C}\text{(}summerseason)\end{array}\right.$ (20)

Figure 7 shows the layout of the CGS plant from which the annual operating data were extracted. The system includes the High Pressure (HP) inlet, two redundant lines with gas Filters (F), Preheaters (PH), gas expansion Valves (V) and the stations for Fiscal Measurement (FM) and Odorant (OD) injection before the gas is fed back into the low pressure (LP).

Figure 7.

City gate station standard layout: red and blue lines represent hot and cold water, dark blue lines natural gas

The control logic of the installed system does not provide for an inverter-controlled flow rate, as assumed in the theoretical study, but a constant value of the pre-heating water flow rate regardless of the gas conditions at the CGS inlet. It affects pre-heating efficiency, as the minimum flow rate of water required for pre-heating is never supplied, and the two pumps (P₁ and P₂) always run at constant revolutions and process the same water flow rate. The boilers operate alternately: the one that does not operate remains on stand-by and consumes an almost constant amount of gas for the pilot flames of approximately 0.25 Sm³/h.

The model used to calculate pre-heating gas consumption was refined by adding improving assumptions. The hypotheses that have been gradually added are as follows:

• Hp₁: Variable gas flow rate crossing the CGS instead of a single constant value,

• Hp₂: Real outlet temperature set point according to DSO,

• Hp₃: Variable gas inlet pressure,

• Hp₄: Variable gas inlet temperature depending on ambient air temperature with the model of equation (9).

Another important hypothesis concerns the control logic of the pre-heating system: to generalise the work as much as possible and after talking with the DSO, it was decided to follow a logic based on modulation of the water flow rate according to the heat supplied to the gas.

Once the water flow temperature to the PH is fixed and the gas conditions are known, the flow rate is derived accordingly to match the thermal demand at PH perfectly.

Figure 8 and Figure 9 present the results of comparing the hourly trends of the thermal load output calculated with the model throughout the year and during four typical days, respectively. The demand curves of the model and the real data are dimensionally adjusted for confidentiality reasons and to make the treatment generic for any CGS: the model with real inputs can faithfully replicate consumption trends for all seasons of the year, as seen in the abovementioned figures.

Figure 8.

Hourly required dimensionless thermal power: model vs. real data during one year

Figure 9.

Daily comparison between the hourly required dimensionless thermal power: final model vs. real data for four different seasons

Figure 9 shows that the model is very accurate during system operation at maximum load (peak daytime hours and winter seasons). On the other hand, the model tends to under- or over-estimate at other times of the day, particularly at night and in summer, when the gas flow rate is very low, and the plant presumably retains a certain amount of heat loss. It can be seen from the figures that the operation of the cabin is highly influenced by the OCs in which it operates: when it is at maximum load, the relationship in equation (1) between consumption W_gas and gas flow rate Q_gas, with all the assumptions added during the development of the model, is optimal.

In the intermediate seasons (autumn and spring), there is a slight underestimation of the typical days shown in Figure 9. For summer, when the gas inlet temperature is very high, the model still estimates a certain theoretical percentage of gas to be burned when it is possible that in real OCs, the cabin exploits thermal inertia to avoid turning on the boilers at times of lower demand. The demand curve used in the following paragraphs is obtained by multiplying the validated dimensionless curve of Figure 8 by the peak power value of 28,738 kW.

Table 2 enlists each improving hypothesis assumed to compute the annual thermal energy consumption with the proposed model. Including a realistic, not estimated, constant flow rate throughout the year, the performance of the thermodynamic model increases significantly (from 413 MWh to about 101 MWh). The following hypotheses allow us to reduce the forecasting from 101 MWh to 74 MWh (Hp₂), then to 53 MWh (Hp₃) and finally to reach a value of E_gas,y of about 44 MWh.

The PV field simple model is taken from the Energy Plus 8.0 database, and thus the every-hour electrical power is given by the following equation:

$W_{PV}(t)=G(t)\times A_{PV}\times {\eta }_{PV}\times {\eta }_{loss}\times f_{cell}$ (21)

Where G(t) is the total solar incident irradiation on the solar panel [W/m²], A_PV is the solar panels’ total area [m²], η_PV and η_loss are the efficiencies of the panels and the inverter system, respectively, and f_cell is the fraction of usable cells. For this work, the three latter parameters are assumed constant and equal to 0.2, 0.98, and 0.95.

The HP is a very efficient technology for heating and cooling purposes since its efficiency, expressed by the COP (Coefficient of Performance), usually varies from 2 to 5 and is particularly high when used to heat a utility or process. The value of the COP is calculated according to (22), considering the efficiency dependence on the temperature difference between the water supply temperature and the ambient air temperature ΔT_HP [22].

$CO{P}_{ref}(t)=6.81-0.121\times \Delta T_{HP}(t)+0.00063\times \Delta T_{HP}{\left(t\right)}^2$ (22)

$\Delta T_{HP}(t)=T_{wat,out}-T_{env}(t)$ (23)

The real COP is then obtained by a correction on the previous formula, which is given for reference value of COP_ref equal to 3.9, including the value of the ideal COP of the considered HP model COP_des.

$COP(t)=\frac{CO{P}_{ref}(t)}{CO{P}_{ref}}\times CO{P}_{des}$ (24)